function [theta1,Q] = LMoptimizerTheta1(theta1Values0,setup,lambda0,tolFun,tolX,maxIter,display)
%This function minimizes the objective function in step 1 with respect to Theta1 using the 
%the Levenberg Marquardt algorhtm. If we find an undefined region then we
%perturb the direction of the search. 

% Locating missing data observations
selectMat       = ~isnan(setup.yields');
totalObs        = sum(sum(selectMat));
theta1          = theta1Values0;
numTheta1       = size(theta1,1);
[Q,output,errorMes] = objectFunc(theta1,setup);
if errorMes ~= 0
    disp('LM optimizer: Initial point is not defined')
    return
end
epsValue        = 1e-8;
residuals       = reshape(output.residuals(selectMat),totalObs,1);
diffQ           = 1e6;
idxIter         = 1;
lambda          = lambda0;
updateJ         = 1;
J               = ones(size(residuals,1),numTheta1);
dThetaMax       = 1;
while (abs(diffQ) > tolFun || dThetaMax > tolX) && idxIter < maxIter 
    if updateJ == 1
        % The gradient
        [~,J] = nlsTheta1(theta1,setup);
        J = -J;
        % Note res = y - yfit. So Grad(res) = -Grad(yFit)
        % residuals_p_eps = zeros(totalObs,numTheta1);
        % errorMesJ = 0;
        %for i=1:numTheta1
        %    theta1_p_eps              = theta1;
        %    theta1_p_eps(i,1)         = theta1_p_eps(i,1) + epsValue;
        %    [~,output_p_eps,errorMes] = objectFunc(theta1_p_eps,setup);
        %    errorMesJ                 = max(errorMes,0);
        %    if errorMesJ == 0 
        %        residuals_p_eps(:,i)  = reshape(output_p_eps.residuals(selectMat),totalObs,1);
        %    end
        %end
        % We have no errors and the gradient is updated
        % Note: Res = y-g(x) so, dRes = -Dg(x)
        %if errorMesJ == 0
        %    J = (-residuals_p_eps - repmat(-residuals,1,numTheta1))/epsValue;
        %end
    end
    
    % Computing the new value of theta1
    %delta = (J'*J + lambda*diag(diag(J'*J)))\J'*(residuals);
    delta = (J'*J + lambda*eye(numTheta1))\J'*(residuals);
    
    % Evaluating fit at new value of theta1
    theta1_new    = theta1 + delta;
    idxSearch     = 0;
    if any(isnan(theta1_new)) || any(isinf(theta1_new)) 
        errorMes   = 1;
        theta1_new = theta1;
        delta      = ones(numTheta1,1);
    else
        [Q_new,output_new,errorMes] = objectFunc(theta1_new,setup);
        diffQ     = Q_new - Q;
        while diffQ > 0 && idxSearch < 30
            idxSearch = idxSearch + 1;
            theta1_new = theta1 + 0.5^idxSearch*delta;
            [Q_new,output_new,errorMes] = objectFunc(theta1_new,setup);
            diffQ  = Q_new - Q;
        end
    end
    perturb = 0;
    indexPer = 0;
    while errorMes ~= 0 && indexPer < 100
        % we perturbe the direction of the search
        indexPer = indexPer + 1;
        perturb = 1;
        theta1_new  = theta1 + delta.*(randn(numTheta1,1).*sum(abs(delta)));
        [Q_new,output_new,errorMes] = objectFunc(theta1_new,setup);
    end
    [dThetaMax,theta1Pos] = max(abs(theta1_new  - theta1));
    if display == 1
        disp(['dQ        = ', num2str(Q_new - Q)]);
        disp(['Q         = ', num2str(Q_new)]);
        disp(['dThetaMax = ', num2str(dThetaMax),' (at element ',num2str(theta1Pos),')']);
        disp(['lambda    = ', num2str(lambda)]);
        disp(['updateJ   = ', num2str(updateJ)]);
        disp(['idxSearch = ', num2str(idxSearch)]);
        disp(['perturb   = ', num2str(perturb)]);
        disp(['iteration = ', num2str(idxIter)]);
        disp(' ');
    end
    if Q_new < Q
        theta1 = theta1_new;
        residuals = reshape(output_new.residuals(selectMat),totalObs,1);
        Q = Q_new;
        updateJ = 1;
    else
        updateJ = 0;
        lambda = lambda*2;
    end
    if idxSearch == 0
        lambda = lambda/2;
    else
        lambda = lambda*2;
    end
    if dThetaMax < 1e-12
        if display == 1
            disp('Change in x is less than 1e-12, so we stop the optimization');
        end
        diffQ = 0;
    end
        
    % Update index
    idxIter  = idxIter + 1;
end
if display == 1
    disp('END OF LM OPTIMIZER')
end
end


